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ABSTRACT 

A recently developed spectral-element adaptive refinement incompressible magnetohydro- 
dynamic (MHD) code [Rosenberg, Fournier, Fischer, Pouquet, J. Comp. Phys. 215, 59-80 
(2006)] is applied to simulate the problem of MHD island coalescence instability (MICI) in 
two dimensions. MICI is a fundamental MHD process that can produce sharp current layers 
and subsequent reconnection and heating in a high-Lundquist number plasma such as the solar 
corona [Ng and Bhattacharjee, Phys. Plasmas, 5, 4028 (1998)]. Due to the formation of thin 
current layers, it is highly desirable to use adaptively or statically refined grids to resolve them, 
and to maintain accuracy at the same time. The output of the spectral-element static adaptive 
refinement simulations axe compared with simulations using a finite difference method on the 
same refinement grids, and both methods are compared to pseudo-spectral simulations with 
uniform grids as baselines. It is shown that with the statically refined grids roughly scaling 
linearly with effective resolution, spectral element runs can maintain accuracy significantly 
higher than that of the finite difference runs, in some cases achieving close to full spectral accuracy. 



1. Introduction 

In many hydrodynamic or magnetohydrody- 
namic (MHD) applications in astrophysics or 
space physics, it is essential that a numerical 
simulation resolve the development of sharp spa- 
tial structures accurately. While pseudo-spectral 
methods generally can maintain high accuracy, 
they are mainly applied on more regular geometry 
and require more uniform grids, which can make 
it difficult to reach high resolution in order to re- 
solve sharp isolated structures especially in flows 
dominated by such structures. Static or adaptive 
mesh refinement (AMR) methods can put more 
grid points in and around isolated structures in 
order to resolve them, but may not achieve similar 
high-accuracy if low order spatial discretizations 
are used ([Rosenberg et al. 2006p. They are par- 



ticularly useful in bounded fiows, where pseudo- 
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spectral methods are often not optimized. There- 
fore, it is of great interest to develop numeri- 
cal schemes that can combine high accuracy and 
high spatial resolution. Adaptive spectral element 
methods have the potential to do just that, pro- 
viding spectral-like accuracy that can be applied 
efficiently to resolve isolated structures. In this 
paper, we concentrate on comparing the accu- 
racy of simulation results from a spectral element 
based AMR code ( SE) ( [Rosenberg et al. 2006| 
[Rosenberg et al. 2006| ) and a finite difference 
based AMR code (FD) ([Germaschewski et al. 20061 
[Bhattacharjee et al. 2005 1 on an astrophysical 
problem that requires high spatial resolution 
as well as high accuracy, the so-called MHD 
island coalescence instability (MICI) problem 
( [Ng fc Bhattacharjee 1998| ). In order to make 
meaningful comparisons, we let each code run on 
essentially the same non-uniform grid (i.e. with 
the total degrees of freedom in the problem fixed) 
that is refined a priori in regions of the grid where 
the current sheets will form in the MICI. A sep- 
arate pseudo-spectral code (PS), running on uni- 



form grids, is also used to provide a baseline for 
the comparisons. 

There is little in the literature about attempts 
to compare FD and high order methods. Often, 
it is simply accepted that the fixed lower order 
FD methods will be less accurate, but that, due 
to their simplicity and efficiency, /i-type grid re- 
finement can always be carried out to improve the 
accuracy for any real problem. We do not attempt 
in this work to compare performance metrics of the 
methods, preferring instead to focus on their abil- 
ity to produce accurate solutions. To an extent, 
then, it is clear that /i-refinement will improve 
solution accuracy for many problems; however, 
previous work ( [Rosenberg et al. "20 06) suggests 
that local high order may be required in certain 
instances. We consider this issue of /i-refinement 
in FD briefiy in this work. Previous comparative 
work in one dimension (jBasdevant et al. 1986^ 
demonstrates that spectral methods, including 
SE, are more accurate than FD methods due 
primarily to dispersion problems in the latter. 
Nevertheless, this work concludes that the SE 
method is not well suited to the calculation 
of thin internal layers (structures), particularly 
when their location is unknown. This assertion 
is made because, while accuracy is found to be 
good, Basdevant et al. 1986 find that polyno- 
mial orders are required to be inordinately high 
even in one dimension. This conclusion is re- 
futed in later work ( [Mavriplis 1994D , which found 
that the SE method is indeed well suited to this 
type of problem provided adaptivity is used (see 
also ([Rosenberg et al. 2006] )). More recent work 
(St-Cyr et al. 20071 also demonstrates that the 



SE method exhibits smaller errors than a low- 
order finite volume scheme at comparable resolu- 
tions for nearly all of a series of tests of a shallow 
water model on a cubed-sphere grid that is either 
adaptively or statically refined. However, to the 
best of our knowledge, ours is the first work that 
performs a quantitative comparison of FD with SE 
methods in divergence-free nonlinear fiows in two 
space dimensions with nonconforming (statically) 
adaptive grids. 

We continue our introductory remarks on the 
MICI in section [2] where we present details of the 
problem, and provide a motivation for our work. 
In the process, we discuss important properties of 
the MICI as well as some aspects of MHD flows 



that will serve our later discussion. We provide de- 
tails about the simulation set up as well as some 
diagnostic measures used in the comparisons in 
section [Sj In section [4] we present the numeri- 
cal methods used in sufficient detail to elucidate 
the results. Our results are presented in section 
[5] where we offer several types of accuracy com- 
parisons relevant to MHD flows, and the MICI, in 
particular. We conclude in section [6] with a sum- 
mary of our flndings, and some discussion about 
the relative advantages of the methods. 

2. MHD Island Coalescence Instability 

In this section, we provide a brief background 
for the MICI problem that motivates us to perform 
our comparative study in this particular simula- 
tion conflguration. 

It is well known that a substantial part of our 
universe is composed of systems of plasmas, ion- 
ized gases, and conducting fluids. Magnetic fields, 
both fluctuating and large scale, play an impor- 
tant role in the physics of these systems. In 
fact, large-scale magnetic fields have long been 
observed to exist in the solar corona (stellar 
coronae, see e.g. (jParker 1979^ and references 
therein), interstellar space (within a galaxy, see 
e.g., (IFor man et al. 19851 IGrimes et al. 2005p ). 
galaxy clusters (see e.g., ( [Kellogg et al. 1973] 
ISarazin 1986^ 1. We consider the representation 
of magnetohydrodynamics (MHD), as a starting 
point of discussion. For simplicity, we consider 
the incompressible MHD equations with constant 
mass density pQ, 

5fU + u • Vu = -Vp-F j X b-h i^V^u, (1) 

a^b = V X (u X b) + TyV^b (2) 

V-u = 0, V-b = (3) 

where u and b are the velocity and magnetic field 
(in Alfven velocity units, b = B/^/^opo with B 
the induction and /xq the permeability) ; j = V x b 
is the current density; p is the pressure divided by 
the mass density; and now the normalized viscos- 
ity v is basically the inverse of Ry, and resistivity 
rj is basically the inverse of S. 

In a dimensionless form, in which all physi- 
cal quantities are measured by their typical val- 
ues, the MHD equations generally have dissipation 
terms involving higher spatial derivatives of the 
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field quantities with strength characterized by the 
inverse of dimensionless parameters: Rv = VL/v 
is the Reynolds number (with V a typical flow 
speed, L a typical length scale), S — VaL/'H is 
the Lundquist number (with Va a typical Alfven 
speed). We can immediately see that for most as- 
tronomical length scales L, both Ry and S are very 
large numbers such that the dissipation terms in 
the MHD equations can be thought usually to be 
ignored, (i. e., the ideal MHD equations), except 
possibly in regions where there exist steep spatial 
gradients. 

In ideal MHD, in which v = rj = in Eqs. Q- 
([3]), several quadratic quantities, namely energy, 
magnetic helicity and cross helicity are conserved 
exactly in three dimensions. Moreover, a magnetic 
field line is carried by the flow velocity so that the 
topology of the magnetic field configuration can- 
not be changed in ideal evolution. Thus, many im- 
portant physical processes, such as magnetic field 
line reconnection, local fiuid heating, and parti- 
cle acceleration due to parallel electric fields, are 
disallowed. In nonideal MHD, such processes can 
occur within boundary layers, which are regions of 
high spatial gradients in the current density and 
vorticity, described by singular perturbation the- 
ory. The tendency for formation of current and 
vortex singularities in the ideal equations, if and 
when it occurs, is a phenomenon of great inter- 
est because the sites of singularity formation are 
precisely the sites where astrophysically significant 
physical processes such as heating and particle ac- 
celeration can occur. 

E. N. Parker has argued for over three decades 
that current sheets, or tangential discontinuities 
of the magnetic field, do generally exist in a 
magnetic equilibrium (jParker 19721 IParker 19791 
IPaxker 1994jl . If Parker is correct, then the for- 
mation of current sheets, and the subsequent fiuid 
heating and particle acceleration from the release 
of magnetic energy due to the rapid dissipation of 
these sheets can have very significant astronomi- 
cal consequences. One observable consequence is 
the production of X-rays. The solar corona was 
the earliest astronomical object that was observed 
to be emitting X-rays. Based on these observa- 
tions, the temperature of the solar corona is esti- 
mated to be of the order of a million degrees with 
peak radiation at wavelengths about 30 A (i. e., in 
the soft X-ray range). While there may be many 



different heating mechanisms involved in heating 
the corona, it is argued that magnetic fields must 
play an important role among these processes. See 
dNg fc Bhattacharjee 2 008) for our recent discus- 
sion on solar coronal heating theory and more ref- 
erences, e.g., (jKlimchuk 2006 p . 

In Parker's model, a solar coronal loop is 
treated as a straight ideal plasma column, bounded 
by two perfectly conducting end-plates represent- 
ing the photosphere. The footpoints of the mag- 
netic field in the photosphere are frozen (line- 
tied). Initially, there is a uniform magnetic field 
along the z direction. To simplify the consider- 
ation, we may keep the footpoints of the mag- 
netic field on one of the plates (z = ) fixed, 
while the footpoints on the other plate (z = L) 
are subjected to slow, random motions that de- 
form the initially uniform magnetic field. The 
footpoint motions are assumed to take place on 
a time scale much longer than the characteristic 
time for Alfven wave propagation between z = 
and z = L, so that the plasma can be assumed to 
be in static equilibrium nearly everywhere, if such 
equilibrium exists, during such random evolution. 
For a given equilibrium, a footpoint mapping can 
be defined by following field lines from one plate 
to the other. Since the plasma is assumed to obey 
the ideal MHD equations, the magnetic field lines 
are frozen in the plasma and cannot be broken 
during the twisting process. Therefore, the foot- 
point mapping must be continuous for smooth 
footpoint motion. (jParker 1972|) claimed that if 
a sequence of random footpoint motions renders 
the mapping sufficiently complicated, there will 
be no smooth equilibrium for the plasma to relax 
to, and tangential discontinuities of the magnetic 
field (or current sheets) must develop. 

Parker's claim has stimulated considerable 
debate ([van Ballegooijen 1985| lAntiochos 19871 



IZweibel fc Li 19871 [Longcope fc Strauss 1994[ [Cowley et al. 1997 1 
that continues to this day. For review and exten- 
sive references, see (jLow 1990( [Browning 1991| 
IParker 1994^ . The first significant objections to 
Parker's claim of non-equilibrium was raised by 
van Ballegooijen (1985), who argued that smooth 
equilibria must always exist as long as the foot- 
point motion is smooth (or continuous). 

van Ballegooijen (1985) developed his argument 
on a reduced form of the MHD equations (referred 
to hereafter as the RMHD equations) originally 



3 



derived by Strauss (1976) for a low f3 plasma. 
These equations, which also provide the basis for 
many later developments in this work, are: 

dtn + [(i3,n] = d,j + [A,j] + i^\7ln, (4) 

dtA+[4,,A]^d,(b + r^WlA, (5) 

where A is the flux function so that the mag- 
netic field is expressed asb = z + W ±A x z, 
is the stream function so that the velocity field 
is expressed as u = Vj_(j) x z, = — V^^c/) is 
the z-component of the vorticity, J ~ —W']_A 
is the z-component of the current density, and 
[</),j4] = dy(j)dxA — dyAdx4>- An ideal magneto- 
static equilibrium solution of Eqs. Q and ([s]) is 
obtained by setting all explicitly time-dependent 
terms, as well as (/) and 77 to zero. We then obtain 

d,J+[A,J]=Q, (6) 

which can also be written as b • V J = 0. This im- 
plies that the current density J must be constant 
along a given magnetic field-line in an ideal static 
equilibrium. 

Based on this set of equations, Longcope and 
Strauss (1994) argued that even when the mag- 
netic equilibrium is unstable, e.g. due to MICI, it 
will only relax to a second equilibrium (assuming 
that it exists) with very thin current layers with 
thickness less than about 10"'' of the large scale. 

However, another point of view was raised by 
Ng and Bhattacharjee (1998), who argued based 
on a mathematical theorem on the RMHD sys- 
tem that for a given fixed footpoint mapping be- 
tween z = and z — L, there exists only one 
smooth equilibrium. This means that an unstable 
equilibrium will relax ideally to a final state with 
current sheets (tangential discontinuities). This 
scenario has very different implications than those 
predicted by the Longcope and Strauss 1994, since 
energy dissipation and other energetic effects can 
be much stronger for the case with current sheets, 
than the case with smooth but thin current layers. 

Therefore, it is very important to determine 
which of these two scenarios should actually occur. 
However, due to the fact that the current layers 
predicted by ( [Longcope fc Strauss 1994| ) are very 
thin, it is beyond our computational ability if the 
full simulation volume is to resolve to the same 
small scale. To provide a resolution to this prob- 
lem with current computer architectures, one may 



need to apply AMR techniques that put more grid 
points in the regions where distinct structures ap- 
pear. At the same time, to ensure that any such 
numerical study is actually representative of the 
true dynamic solution, one needs to make sure 
that the numerical scheme used can maintain a 
reasonably high accuracy. Hence, accuracy be- 
comes a particularly important factor in the choice 
of the numerical method for the Parker problem 
of MICI. 

Because finite difference-based schemes are 
usually of low order truncation, accuracy typ- 
ically decreases as higher order derivatives are 
taken. Also because of low order, these methods 
can be diffusive (and dispersive). For the RMHD 
equations Q and ([S]) , one needs to use the spatial 
derivatives of J, or third order spatial derivatives 
of A. A couple of questions arise in simulations 
based on such schemes: Will the lower accuracy in 
calculating these higher spatial derivatives change 
drastically the dynamical properties of the prob- 
lem? Will the numerical diffusion preclude the 
formation of a current sheet, and lead instead to 
a smooth residual current layer? For example, in 
the Parker problem, one will need a simulation 
that is accurate enough so that we can show con- 
fidently whether there is a formation of a true 
current sheet as predicted by Ng & Bhattachar- 
jee (1998), or if the current layers actually tend to 
fixed (but small) thickness (as small as 10"'' of the 
large scale) as predicted by Longcope & Strauss 
(1994). In order to increase the reliability of the 
simulations, it is thus of great interest to look for 
schemes that can maintain higher accuracy, and 
can make use of AMR techniques to resolve sharp 
features at the same time. Spectral element meth- 
ods have the potential to fulfill such requirements. 
This is what motivates us to perform the present 
comparative study. 

3. Simulation set up and diagnostics 

While the main problem of interest is the three 
dimensional (3D) MICI problem represented by 
Eqs. (|4])-([6]), for simplicity we will instead simu- 
late the two dimensional (2D) version of the prob- 
lem for the purpose of this comparative study. 
Examining this problem in 2D should not affect 
greatly differences in accuracy among different 
codes. The 2D MICI problem can be viewed as 
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the limit of the 3D problem with L ~f oo. In 
this limit, it has been shown that current singu- 
larities must form for the ideal equations {-q = 0) 



( Longcope fc Strauss 1993 1. This is beneficial for 
our present study since we know that there must 
be specific structures, and we know where they 
will appear. We can then compare how well these 
sharp structures are resolved in different schemes. 
Therefore, we will simulate the set of equations 
without the z-dependence, or 

dtVL + [0, n] = [A, J] + (7) 

dtA + [4>, A] = riVlA, (8) 

where the definitions for the variables are the same 
as in Eqs. (|4])-(|6]). 

Periodic boundary conditions are used in both 
the x~ and y— directions. To be specific, the 
simulation domain is set to be < a; < 1 and 
< y < 1. The initial equilibrium is chosen to be 

A{x, y,t = 0) = Aq sin(27ra;) sm{2TTy) , (9) 

with a small initial flow, which in terms of the 
stream function is 

(l){x,y,t = 0) = <j)o[cos{2nx) ~ cos{2TTy)] , (10) 

where Aq — OA and 00 = 0.002. 0o is chosen small 
enough so that there is a clear linear phase in the 
growth of the MICI. Because SE evolves a different 



form for the equations (section 4.1), these initial 
conditions are converted into conditions on b by 
using b = V X Az, and on u by using the relation 
u = V X 0z. 

Fig. [l] shows the contour plots of A and (p at 
t = as given by Eqs. ^ and (10 1. The velocity 



deduced from this stream function, has the initial 
tendency to push the two islands at the lower right 
and upper left toward each other and towards the 
center. Eventually, these two islands will merge 
with each other; hence, the nomenclature island 
coalescence. 

The forms in ^ and ( [To| should preserve addi- 
tional symmetries with respect to the lines x ~ y 
and X = l-y: A{x,y) = A{~x,~y) = A{y,x), 
(l^{x,y) = 0(-x, -y) = -<j){y,x) . Specifically, 
(j) — on the lines x = y and x = 1 — y. These 
symmetries are not incorporated into the simula- 
tion schemes; however, they can provide informa- 
tion about how well a numerical scheme preserves 
them. 



The grids are generated separately for each Rv 
we consider, and shown in the following sections. 
The equivalent resolution 7?.cq for the FD and 
SE simulations is defined as that which would be 
achieved if the most finely resolved subdomain 
covered the entire domain. If Eq is the linear 
number of elements then for SE, equivalent lin- 
ear resolution, TZcq, is computed from the initial 
Eq X Eq grid, the number of refinement levels, 
£, and the expansion polynomial degree, N, such 
that 7^oq = [2'^NEq]. For all SE runs, Eq = 8, 

= 8, and £ varies with R^. For the FD method, 
in order to make more direct comparisons, we use 
the same grids as used in SE, with a 8 x 8 uniform 
grid within each element when 8*^^ order polynomi- 
als are used in SE. For the PS method, a uniform 
grid of TZcq x TZcq collocation points is used. Table 
[l] contains a list of the viscosity (resistivity) used 
and the corresponding 7?,oq. 

In the absence of external forcing, viscosity and 
magnetic resistivity, the 2D MHD equations ^ 
and m have three ideal invariants: the energy 



Ek + E 



M ■ 



(11) 



composed of the kinetic and magnetic energy, the 
L2 norm of the magnetic potential 



M = {A^) , 
and the cross helicity 

H=(uh) 



(12) 



(13) 



For all functions (p in these definitions, we have 

(0)^/</>dx2. 

In 2D assuming bi-periodicity, it is easy to show 
using Eqs. - ([s]) that 

^ = - ^( J2) , (14) 



Table 1: Parameters used in the simulations de- 
scribed in the following sections; TZcq is the linear 
grid resolution, and v (77) is the viscosity (resistiv- 

ity). 





V = Tj 


128 
256 
512 


2 X 10^^ 
1 X 10-3 

3 X lO-'' 
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dM 



and 



dH 

~dt 



(15) 
(16) 



into the d-dimensional variational form of equa- 



Equations (14l-(16l must hold for any 2D pe- 



riodic MHD flow, and they serve as critical diag- 
nostics for any numerical treatment of the incom- 
pressible MHD equations. 

In addition to the conservation laws, there are 
several other quantities of interest in diagnosing 
numerical solutions to MICI as defined later when 
we present our numerical results. 

4. Numerical methods 

We provide in this section the numerical meth- 
ods we use for carrying out the simulations on 
the MICI problem. These methods have been de- 
scribed in detail in other sources, but we provide 
enough explanation so that the simulation setup, 
results, and discussion will be comprehensible, and 
the paper reasonably self-contained. 

4.1. Spectral element method 

This method evolves Eqs. ([l])-(|3]) in time, as 
part of the Geophysics/ Astrophysics Spectral Ele- 
ment Adaptive Refinement (GASpAR) code, and 



has been described in detail in ( IRosenberg et al. 2006^ 
Here we present aspects of that description that 
relate to the results discussed in this work. 
Eqviatio ns ([l|)- (|3|) are solved in Elsasser form 
(|ElsasserT950| : 

9tZ=^-HZ=F-VZ=t+Vp-j/±V2Z±-i/=Fv2z=F = 

(17) 

V-Z± = 0, (18) 

where = u ± b and — ± rj). Thus, 
we solve for u and b in terms of Z*. All com- 
ponents of Ti"^ and those of the velocity and 
magnetic field reside on Gauss-Lobatto-Legendre 
(GL) nodes, while pressure resides on Gauss- 
Legendre (G) nodes. These choices follow from 
the finite dimensional subspaces to which these 
quantities are restricted: Z* (and u and b) are 
expanded in terms of N^^ order GL polynomi- 
als, while p is expanded in terms of — 2'^ or- 
der G polynomials. Substituting these expansions 



tions (17)-(18l on a domain D, and using appro- 



priate quadrature rules, we arrive at a set of semi- 
discrete equations written in terms of spectral el- 
ement operators: 



± 



dZ, 
dt 



D^Z± 



MC^Z± 

^i^±LZ± - 



.^LZJ 



0, 



(19) 
(20) 



for the j component. The variables Z repre- 
sent values of the Z* collocated at the GL node 
points, and are values of the pressures at 
the G node points. Note that, even though the 
continuous equations contain only a single pres- 
sure, Eq. (19) contains one for each Elsasser vec- 



tor because the constraints ( 20 1 are enforced sep- 
arately on them. The operators M, L, and C, 
are the well-known mass matrix, weak Lapla- 
cian and advection operators, respectively {e.g. , 
( [Rosenberg et al. 2006| ), and references therein), 
and Dj represent the Stokes operators, in which 
the GL basis function and its derivative opera- 
tor are interpolated to the G node points, and 
multiplied by the G quadrature weights. All 
d-dimensional operators are computed as tensor 
products of their component ID operators. Note 
that because different expansion bases are used 
for the vector quantities as are used for pressure, 
a staggered grid is implied. Hence, this method 
is referred to as the P^r — PAf-2 formulation. It 
was chosen to prevent spurious pressure modes 
dMaday et al. 1992[ IFischer 1997)l . Note also the 



effect of the Dj , which act on so-called v-grid (vec- 
tor) quantities: they take a derivative that itself 
resides on the G nodes; hence, the discrete diver- 



gence, Eq. (20 1 resides on the p-grid. The effect of 



the transposed Stokes operator, , on the other 
hand, is to compute a derivative-residing on the 
v-grid-of a p-grid quantity. 

The code has an adaptive mesh capability that 
we utilize minimally in this work. The connec- 
tivity and adaptivity algorithms were described in 
( Rosenberg et al. 2006| ). No dynamic adaptivity 
is used here (see section [s]). Instead, the noncon- 
forming grid is constructed initially by turning off 
the refinement criteria, and selecting the elements 
we want to refine explicitly. The nonconforming 
grid is then used in a static configuration through- 
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out the simulation. Nonconforming in this context 
means that there are at most two child elements 
adjacent to a coarser neighbor; an element is re- 
fined by dividing it isotropically into 2x2 child 
elements, each of which contains the same polyno- 
mial order as the parent. 

4.. 1.1. Time stepping with a constraint 

In our simulations we must resolve all tem- 
poral (and spatial) scales, so we use an ex- 
plicit Runge-Kuttta (RK) method for integrating 
Eqs. p9l)-(p0| in time. The specific RK algo- 



rithm is that presented in (jCanuto et al. 19881 p. 
109), known to be valid to second order in At 
(|Brachet et al. 2007p . which we use for all com- 
putations. At each stage, we can write (recall eq. 



(191): 



(21) 



i/±LZ± + i.^LZj). 



We require that each RK stage obey eq. ( pO| in 
its discrete form, so multiplying eq. (21 1 by Dj, 



summing, and setting the term D-'Zj = leads 
to the following pseudo-Poisson equation for the 
pressures, p"*": 



(22) 



where the quantity 



= l^i M-i(MC^Z±+.±LZ±+.^LZj)-Zf " 

is the remaining inhomogeneous contribution. 
Rosenberg et al. (2007) describe how is 
computed. Equation ( p2| ) is solved using a pre- 
conditioned conjugate gradient method (PCG); 
for all computations reported here, we use a block 
Jacobi preconditioner computed using a fast diag- 
onalization method. 

Thus, at each time-step, two RK stages are 
computed, and each stage solves eq. (22) twice, 
once for Z+, and once for Z~ leading to 4 pseudo- 
Poisson solves at each time step. 

4.1.2. Communication 



While equations (19l-(20l are correct for a 
single element, strictly speaking, they are not 
complete when multiple subdomains (elements) 



are used. In this case, we must ensure that all 
quantities in the subspace represented by the GL 
grid remain continuous across element interfaces. 
The specific way in which this is carried out in 
the code is described in ( [Rosenberg et al. 2006| ), 
and entails exchanging interface data. Using the 
Boolean scatter matrix, Ac, the interpolation ma- 
trix from global to local degrees of freedom, 
and the masking matrix (that enforces homoge- 
neous boundary conditions), fl, that were pre- 
sented there (see also Fischer et al. 2002 ), it is 
found that the local Stokes operators in equations 



(19l-(20l must be replaced with 



Dj -> Dij*Aen, 

where D^j — diag;.(D^), and the are the ma- 
trices from above, and k indexes the elements in 
the domain. It is sufficient for our purposes to re- 
fer to the local form of the Stokes operators, and to 
simply observe that communication occurs when 
they are applied. 

4.2. Finite difference method 

The Magnetic Reconnection Code (MRC) is 
a suite of codes (jCermaschewski et al. 2006] 
[Bhattacharjee et al. 2005| which integrate vari- 
ous reduced and extended fluid models of plasma 
flows. In this paper, the 2D AMR version of the 
code integrating the equations of RMHD has been 
used. 

The MRC employs a hierarchical quadtree 
based approach in block-structured adaptive 
mesh reflnement. At each reflnement time, ev- 
ery (square or rectangular) box is checked as to 
whether the data in that box requires additional 
resolution. The reflnement criterion needs to be 
selected appropriately for the given problem, from 
simple evaluation of maxima or gradients to a 
Richardson extrapolation based approach. 

If the local resolution in a box is considered in- 
sufficient, it is subdivided into 2x2 smaller boxes, 
which have physically twice the resolution but the 
number of grid points in each box remains the 
same it was on the coarse parent box. 

As opposed to the alternate approach of allow- 
ing for arbitrary rectangular patches of reflnement, 
this method results in simpler data management, 
easier optimization (each box is always the same 
size of n X m grid points, important for cache con- 
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siderations, etc) and more efficient load balancing, 
even though it introduces some inefficiency in that 
some additional areas are refined where the higher 
resolution is not required. 

For load-balancing, we use a space-filling 
Hilbert Peano curve which connects all the boxes 
at various levels of refinement. The boxes along 
this one-dimensional curve are then evenly dis- 
tributed to the available processors. Since each 
patch has the same number of grid points, the 
computational load is evenly spread, and as the 
Hilbert-Peano curve has the property that in a 
certain averaged sense, patches which are close 
in the two-dimensional domain are also close on 
the one-dimensional curve, spatially close regions 
are clustered onto the same processor; that is, it 
maintains data locality. Since communication is 
only necessary between spatial neighbors to ex- 
change boundary data, most of the needed data 
is available on the local processor, expensive MPI 
communication is only required for boundary data 
transfer across the boundaries between clusters 
on different processors, which is effectively mini- 
mized. 

An additional difficulty in using AMR to in- 
tegrate the reduced models - not encountered 
in the purely hyperbolic systems - is the need 
to solve elliptic equations, e.g. solving for the 
stream function (j) from the vorticity f2. To dis- 
cretize this sub-problem, we rewrite the Lapla- 
cian as the divergence of a gradient V • V 
and apply a conventional finite volume method 
for a conservation law. Combined with correcting 
the fluxes at fine-coarse boundaries, this provides 
an exactly numerically conservative discretized ex- 
pression. To efficiently solve this elliptic prob- 
lem we use a variation of the fast adaptive com- 
posite (FAC) method (jGermaschewski et al. 20041 
IMcCormick 1989)1 . which, due to the multilevel 
character, provides an iterative solver with very 
fast convergence. 

Alternatively, the discretized Laplacian on the 
AMR hierarchy can be calculated explicitly as a 
sparse matrix, and PETSc's fBalay e t al. 2004^ 
rich supply of solvers are available to solve 
Poisson's equation. In particular, SuperLU 
(|Li k Demmel 2003^ in many cases proves to be 
very fast, since the expensive LU decomposition 
only needs to be done once and then can be ap- 
plied for many time-steps until the AMR hierarchy 



of boxes changes. This is the method used in the 
simulations for this study. 

Our AMR code can also be run in fully implicit 
mode using Crank-Nicholson time-stepping, how- 
ever no preconditioner has been developed yet, so 
we are using a simple unpreconditioned Newton- 
Krylov solver which does not achieve optimal per- 
formance, but is still faster than a Runge-Kutta 
explicit scheme. Most of the runs presented here 
axe using the implicit time integrator. Some runs 
have also been cross-checked by using an adaptive 
time-step Runge-Kutta explicit method, which is 
included with PETSc's. 

Spatial discretization is using a regular second 
order central differences, which is equivalent to a 
conservative finite volume method, with no up- 
winding employed. 

For the purpose of the present study, the dy- 
namic refinement capability of the AMR code is 
turned off so that we can use the same grids that 
were also used in SE. This was done because it is 
difficult to set refinement methods and criteria to 
be the same in both the FD and SE codes. Using 
static refinement provides us with nonconforming 
grids over which we can exert complete control of 
the number of degrees of freedom. Comparison of 
refinement criteria in the SE and FD codes and the 
effect on the resulting dynamics is left for future 
work. 

4.3. Pseudo-spectral method 

The PS code is based on fast Fourier trans- 
form (FFT) on a 2D bi-periodic domain. It is 
de-aliased by the standard 2/3 rule. The nonlin- 
ear term is calculated in the physical space on a 
uniform grid of collocation points. A second order 
predictor-corrector method is used for time inte- 
gration. The code is parallelized using a parallel 
version of the FFT. For this study, results from SE 
and FD methods for the TZcq case are compared 
with those from the PS code with TZ^q x TZ^q col- 
location points. The results from the PS code are 
themselves checked with PS runs with higher res- 
olutions, up to 2048 X 2048, which confirm that 
the results with the original resolution are already 
well resolved. 
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5. Numerical results 




0,0 0,2 0,4 J O.L) 0.(1 1.0 0.0 O.a 0.4 3; 0.6 O.S 1.0 

(a) (b) 



Fig. 1. — Contour plots of (a) the initial flux func- 
tion A, and (b) the initial stream function 0. Pos- 
itive (including zero) contours are solid and nega- 
tive contours are broken. Contours levels are from 
—0.4 to 0.4 with an increment of 0.025 in (a), and 
-0.004 to 0.004 with an increment of 0.00025 in 
(b). 



A{1 = 1) 0(1 = 1) 




(«) (b) 



Fig. 2. — Contour plots of (a) A, and (b) at 
t = \ for the 7?.eq — 256 case. Positive (including 
zero) contours are solid and negative contours are 
broken. Contours levels are from -0.4 to 0.4 with 
an increment of 0.025 in (a), and -0.1 to 0.1 with 
an increment of 0.00625 in (b). Note the sheet 
formation (a) and the fluid acceleration (b). 



In this section, we compare simulation solutions 
by SE and FD methods, using solutions by PS 
as references. Note that each code (SE and FD) 
has been verified separately on a suite or problems 
including analytical solutions, but the purpose of 
the present work is to explore the fully developed 
nonlinear regime of the MICI problem, when sharp 
current and vorticity layers appear. 

The main dynamics of MICI is due to the at- 
tractive force between two islands with the same 
sign of current. The initial small perturbation 
in the flow velocity breaks the unstable equilib- 
rium. After a short transient, a linear phase ap- 
pears when kinetic energy increases exponentially 
so that the flow velocity pushes the two islands 
together further. Eventually a sharp current layer 
appears between the two islands when the dynam- 
ics enters the nonlinear phase at time t ~ 1, in 
units of the Alfven time. Magnetic reconnection 
then proceeds faster, producing a strong outflow 
velocity, as well as sharp vorticity layers. Fig. [2] 
shows contour plots of A and (/> at f = 1 for the 
7?.eq = 256 case. We see that the two islands with 
positive A are pushed towards each other with 
some A contours (magnetic field lines) already re- 
connected, as compared to Fig. [T] The stream 
function shows strong outflows as indicated by 
the concentration of stream lines. 

The grids used in our simulations for the three 
cases are shown in Fig. [3] to |5] The reflnements 
with each grid are imposed so as to resolve struc- 
tures produced by the above dynamics, for the 
different dissipation levels indicated in Table [l] 
In Fig. |3] we also plot the color contours of the 
current density J at i = 1.3 produced by the SE 
run for the 7?,oq = 128 case, with red representing 
the positive end and blue representing the nega- 
tive end of J values (as in the following). Only 
one level of reflnement is used that covers the re- 
gion containing stronger J. Note that within each 
square, polynomials of order = 8 are used for 
SE, while an 8 x 8 uniform grid is used for FD 
(the same as in the other two cases). In Fig. |4] 
the color contours of the vorticity at i = 1.3 
produced by the SE run are plotted along with the 
grid for the 7?.oq = 256 case. Note that there are 
now two levels of refinement. In Fig. |5] the color 
contours of J at f = 0.93 produced by the SE run 
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Fig. 3. — The statically refined grid used in the 
TZcq = 128 case. The background is the color con- 
tour plot of the current density, J, at t — 1.3 pro- 
duced by the SE run, with red representing the 
positive end and blue representing the negative 
end of J values. Within each square, polynomials 
of the order of iV = 8 are used for SE, while an 
8x8 uniform grid is used for FD. Only one level of 
refinement is needed at this (low) Reynolds num- 
ber. 
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Fig. 4. — The grid used in the 7?.oq = 256 case. 
The background is the color contour plot of the 
vorticity at t = 1.3 produced by the SE run, 
with red representing the positive end and blue 
representing the negative end of values. Two 
levels of refinement are used here. 
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Fig. 5. — The grid used in the 7?.cq — 512 case. 
The background is the color contour plot of the 
current density J at t = 0.93 produced by the SE 
run, with red representing the positive end and 
blue representing the negative end of J values. 
Now, three levels of refinement are used. 
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are plotted along with the grid for the TZ^q = 512 
case, with one further level of refinement. At the 
same time, we found that the largest squares have 
to be refined for this case due to a much smaller 
dissipation in this case. 

Before presenting our numerical results, let us 
look at the degrees of freedom (DOF) of these 
grids at the diS'erent TZeq levels, as shown in Fig. [6] 
in red asterisks; they follow roughly a linear scal- 
ing proportional to TZcq- Also plotted are DOF of 
the PS runs using uniform grids, in blue diamonds, 
which follow a scahng of TZlq as they should. Since 
as TZeq increases, the diS'erence between the two 
scalings can be very large, using adaptive grid re- 
finement has the potential to provide considerable 
savings in memory and/or CPU usage, if the linear 
scaling of the DOF in these adaptive grids contin- 
ues to hold for even larger TZeq, and thus for large 
Rv and S. 
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5.1. Accuracy of the conservation laws 

We start our comparison by looking at how 
each code preserves conservation laws as shown 



in Eqs. ( 14 1 and ( 15 1. We do not include the con- 



servation law for cross helicity, Eq. (16 1, in our 



comparison, since H — exactly based on our 
initial conditions in Eqs. (l9| and (10 1, and thus 



dH/dt — 0; all three codes preserve H = and 
dH/dt — well during the duration of the runs 
for all three cases. This is more due to how well 
the structure of each code preserves symmetries, 
rather than due to numerical accuracy. Also, since 
both the left hand side (LHS) and right hand side 



(RHS) of Eq. ( 16 ) are close to zero, taking the dif- 



ference between the two to see fractional changes 
will not yield meaningful results. 

In Figs.|7]to|9] the fractional difference between 
the LHS and RHS of (a) Eq. (fTl), and (b) Eq. (Il5l 



are plotted as functions of time, for the three val- 
ues of 7?,eq. The fractional difference (or error) A is 
defined as | LHS - RHS | / | RHS |, with the time 
derivative in the LHS calculated by taking finite 
difference of the output in time, thus providing an 
overestimate of the error. In these figures, black 
curves are for PS runs, red curves are for SE runs 
and blue curves are for FD runs. 

For the TZeq — 512 case, we present results up 
to t ~ 1, since both SE and FD experience larger 
error after this time, probably due to the fact 



Fig. 6. — DOF of the grids shown in Fig. [3] to [5] at 
the different TZeq levels, in red asterisks, and those 
of the PS runs using uniform grids, in blue dia- 
monds. A dashed line showing a TZlq dependence 
in the PS runs, and a dotted line showing a linear 
dependence of TZeq of the statically-refined runs 
are also plotted. 




Fig. 7. — Fractional error in energy conservation 
law in (a) AE, and for the C2 norm of the mag- 
netic potential, AM in (b) for the TZeq = 128 case 
as functions of time. Black curves are for the PS 
run, red curves are for the SE run, and blue curves 
are for the FD run (same in the other figures be- 
low) . 
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(a) (b) 



Fig. 8. — Fractional error in conservation laws in 
(a) AE, and (b) AM for the 7^cq = 256 case as 
functions of time. See Fig. [7] for definitions. 




0,0 0,2 0,4 0,6 0,8 1,0 OS 4 6 OR 1,0 

(a) (b) 



Fig. 9. — Fractional error in conservation laws in 
(a) AE, and (b) AM for the 7^oq = 512 case as 
functions of time. See Fig. [7] for definitions. 



that we are using a fixed adaptive grid and could 
not follow the development of small scales closely 
enough. To achieve more accurate results, both 
codes would have to employ dynamic adaptive re- 
finement. However, it is difficult to make certain 
the two codes refine and coarsen in the same way 
in order to make meaningful comparisons; hence, 
this investigation is left for future work. 

For all three levels of 7?.cq, we see that PS results 
(black) preserve conservation laws the best, as ex- 
pected. This is why we may use them as baselines 
for comparisons. The error level in AM as shown 
in (b) panels (with M = d{A'^)/dt) is around or 
slightly over 10~^, while the error AE as shown in 
(a) panels is around or slightly below 10^^, since 
energy involves one more spatial derivative than 
A and thus is less accurate in this computation. 

For SE runs (red) , the accuracy turns out to be 
quite good, given that they are running on stati- 
cally refined grids. The error AE is more or less 
one order of magnitude above that of PS runs for 
all three cases, but still is in a low level of about 
10"''. For the error AM, the difference between 
SE and PS becomes larger (about two orders of 
magnitude), at a level of about 10~'^. The reason 
AM is greater than AE for SE runs is due to the 
fact that b and u are the primary variables in the 
computations, and A is a field quantity derived 
from b, by solving the equation V^A = — J. This 
process introduces error in A and thus the error 
AM is greater than AE. 

For FD runs (blue), the code does capture 
quantitatively the main dynamical evolution of the 
MICI problem, although with a lower accuracy. 
For AM, the error level is up to above 10~^, and 
the error AE can go above that, even reaching as 
high as 10^^ or more. The error in AM is at a 
level slightly lower than that of in this version 
of the FD code, similar to the trend of PS, since 
it uses A as primary variable. 

From these comparisons of the conservation 
laws, we see that the accuracy of SE is higher 
as expected for spectral-based methods; it some- 
times approaches that of PS runs using uniform 
grids. This confirms that SE can deliver near spec- 
tral accuracy, the main advantage of employing 
such a scheme. Using the same adaptive grid, FD 
runs show lower accuracy. It is conceivable that 
the accuracy of the FD scheme can be improved 
somewhat by algorithm modifications. However, it 
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is unlikely that it can be improved to the level of 
the SE runs, when constrained to the same grids. 

5.2. Accuracy with respect to higher order 
derivatives 

In the above comparison of conservation laws, 
we see that the accuracy level can change with re- 
spect to field quantities involving a different order 
of spatial derivatives of A. We now look into this 
issue further by comparing integrated field quan- 
tities. Instead of presenting comparisons for all 
three levels of 7?.oq, we will only present figures 
based on the TZcq = 128 case, and just mention 
that similar conclusions can be drawn for the other 
two cases. 
In Figs. 

(b2 



;o 13 we show in the (a) panels the 
time series of the quantities (A^), E = ^(b^ + u^), 
(J2), and ((VJ)2), for the PS run (black), SE 
run (red), and FD run (blue), plotted in this or- 
der for the TZcq = 128 case. These four plots in- 
volve field quantities of increasing order of spatial 
derivatives of the magnetic potential A. When 
plotting this way, it is not easy to see the differ- 
ence between the curves when they are close to 
each other. In fact, red curves almost totally cover 
black curves for all cases, and blue curves almost 
totally cover the other two in {A"^) and E. So in 
the (b) panels, we plot the fractional difference 
A between these values and those from a "con- 
verged" PS run using a much higher resolution 
with a uniform grid of 2048 x 2048. Again, black 
is for the PS run (fractional difference between the 
Tleq = 128 PS run and the 2048 PS run), red is 
for the SE run, and blue is for the FD run. The 
fractional difference (or error) A is defined as, e.g., 

A{A^) ^1 (A^) - (^2)2048 I / I {A%04S |. 

We note from the black curves that the PS run 
with a 128 x 128 grid is indeed a highly converged 
run, in the sense that the fractional error with 
respect to the 2048 x 2048 is very small. The er- 
ror A{A'^) is at a level of about or slightly over 
10~^. This increases to around 5 x 10"'' for AE, 
around 10^^ to 10^^ for A(J^), and with a maxi- 
mum value shghtly below 10"^ for A((VJ)^). The 
trend of obtaining less accurate results for quanti- 
ties involving higher derivatives is expected. Still, 
for all quantities that are important to the numer- 
ical integration, i.e. up to V J, the accuracy of the 
128 x 128 run is high enough for the whole dura- 
tion of the simulation. Again, this shows that the 



(a) 




' (b) ' 



Fig. 10.— (a) Time series of {A'^) for the PS 
(black), SE (red), and FD (blue) runs plotted in 



this order for the TZ, 



eq 



128 case. All three runs 



overlap, (b) Fractional error in A(A ), as com- 
pared with a PS run with 2048 x 2048 resolution. 
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t 



Fig. 11. — (a) Time series of total energy E for the 
PS (black), SE (red), and FD (blue) runs plotted 
in this order for the TZeq = 128 case. Again, all 
three runs are nearly coincident, (b) Fractional 
error in AE, as compared with a PS run with 
2048 X 2048 resolution. 




Fig. 12.— (a) Time series of (J^) for the PS 
(black), SE (red), and FD (blue) runs plotted in 
this order for the TZ^q — 128 case. Differences in 
the dissipation ?]( J^) are observed, (b) Fractional 
error in A(J^), as compared with a PS run with 
2048 X 2048 resolution. 
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PS run can be used as a baseline for comparison. 

For the SE run, the error level is higher than 
that of the PS run. This is qualitatively similar to 
the comparison of conservation laws, except now 
it is higher by a somewhat larger amount, some- 
times about two orders of magnitude. Specifically, 
A(j4^) is at a level of about or a little over 10~^, 
is up to about 5 x 10"'', A(J^) is about or a 
little over 10^^, and A((VJ)^) is about or a lit- 
tle over 10~^ . This level is still acceptably low, 
but is somewhat higher than those in the error of 
conservation laws as shown above, especially for 
A(J^) and A((VJ)^). However, we need to keep 
in mind that the errors in the conservation laws are 
indicating how well a code simulates the equation 
self-consistently at each moment, but the errors in 
comparing with results from a converged solution 
are accumulated over time, and can thus be larger 
than the former. 

For the FD run, the error level is again higher 
than that of the SE run, sometimes by an order 
of magnitude or more. The error level of A(A^) 
is about 10""^ or slightly above, while Ai5 is at 
a level or slightly over 10^^. These two are still 
reasonably small and so we do not observe ap- 
preciable differences in their respective (a)-panel 
plots. However, A(J^) and A((VJ)^) become 
large enough to be observable in the time series 
plots themselves. Quantitatively, A(J^) is at a 
level of 10~^ or a slightly below, while A((VJ)^) 
is at a level of 10~^ or somewhat above. Although 
such high error levels do not seem to alter the over- 
all dynamics of the solutions qualitatively, they are 
at a level high enough to be of some concern. 

5.3. Accuracy of current layer width 

One quantity of great interest in the MICI prob- 
lem is the width of current layers, which can be 
defined as 



l^{{J^)/{{VJf)) 



(23) 



In Fig. |14[ we show the current layer width I 



defined in Eq. ^ for the PS (black), SE (red), 
and FD (blue) runs, plotted in this order for the 
TZcci = 128 case in panel (a), and the error AZ, 
as compared with the 2048 x 2048 PS run in (b). 
We see that the error level A/ is similar to that of 
A((VJ)^), at about lO^'^ or slightly over for the 
FD run, about an order of magnitude higher than 



that of the SE run. Again, this accuracy level 
is qualitatively reasonable. However, if we need 
to investigate the more difficult problem of the 
Parker's model (3D MICI) and need to determine 
whether Z ^ or not in the 77 ^ limit, a 10% 
error could lead to significant uncertainty. 

5.4. Quantitative summary 

Let us summarize the comparisons from Fig. [?]- 
|13| by plotting the maximum fractional error over 
the plotted time period for all cases. In Fig. [Ts] 
we show in (a) the maximum fractional error over 
the plotted period in Fig.[7]to Fig.|9]of Ai; (black 
triangles for PS, red asterisks for SE, blue dia- 
monds for FD) and AM (black plus signs for PS, 
red squares for SE, blue crosses for FD) for the 
three 7?,oq cases. In (b), we show the maximum 
fractional error over the plotted period in Fig. [TO] 
to Fig. [l3]of A(A2) (n = 0), A£; {n = 1), A(j2) 
(n = 2), and A((V J)^) (n = 3) versus n, the order 
of spatial derivative with respect to A (black tri- 
angles for PS, red asterisks for SE, blue diamonds 
for FD), for the 7?,eq = 128 case. From these two 
plots, we see clearly the separation between the 
three groups (black for PS, red for SE, and blue 
for FD). While PS has the lowest level of error as 
expected (due to a fully spectral treatment and to 
the fact that uniform grids are used), SE achieves 
an error level somewhat in between PS and FD. 
In some cases, the error in SE approaches that of 
PS, indicating near spectral accuracy. 

So far we have looked at comparisons of spa- 
tially integrated quantities. It is conceivable that 
the accuracy of field quantities at particular spa- 
tial points may not follow similar trends. Here we 
present one example that compares values of the 
maximum current density, Jmaxi over the entire 
box. In Fig. |16| we show the maximum current 
density value over the periodic box, Jmax, for the 
PS (black), SE (red), and FD (blue) runs, plot- 
ted in this order for the TZ^q — 128 case in panel 
(a), and the error AJ^axi as compared with the 
2048 X 2048 PS run in (b). 

We see in this comparison that even the PS run 
has an error level at around 10~^, much higher 
than the errors of other quantities we have shown 
so far with this method. This is because the grid 
used in the 2048 x 2048 run is much finer than 



what is used in the TZ„ 



128 PS run. So, the 



value of Jniax from a much higher resolution run 
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Fig. 13.— (a) Time series of ((V,/)^) for the PS 
(black), SE (red), and FD (blue) runs plotted in 
this order for the TZcq = 128 case. PS and SE 
data still overlap, but FD the FD result is notice- 
ably different, (b) Fractional error in A((VJ)^), 
as compared with a PS run with 2048 x 2048 res- 
olution. 



Fig. 15. — (a) Maximum fractional error over the 
plotted period in Fig. [t] to Fig. [o] of AE (black 
triangles for PS, red asterisks for SE, blue dia- 
monds for FD) and AM (black plus signs for PS, 
red squares for SE, blue crosses for FD) for the 
three T^-oq cases, when compared to a PS run at 
a grid resolution of 2048a; x 2048. (b) Maximum 
fractional error over the plotted period in Fig. [TO] 
~of A(A2) (n 0), AE (n = 1), A(j2) 
2), and A((V J)^) (n — 3) versus n, the order 



to Fig. 

in. 

of spatial derivative with respect to A (black tri- 
angles for PS, red asterisks for SE, blue diamonds 
for FD), for the 7^eq = 128 case. 




Fig. 14. — (a) Time series of the current layer 
width I for the PS (black), SE (red), and FD (blue) 
runs plotted in this order for the TZcq — 128 case, 
(b) Fractional error in Al, as compared with a PS 
run with 2048 x 2048 resolution. 




1 (a) ^ = " ' (b) ^ 



Fig. 16. — (a) Time series of Jmax for the PS 
(black), SE (red), and FD (blue) runs plotted in 
this order for the TZcq — 128 case, (b) Fractional 
error in A Jmax, as compared with a PS run with 
2048 X 2048 resolution. 
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cannot be located at the collocation points of the 
lower resolution run, and thus it reaches a slightly 
higher value. With this consideration, the fact 
that both the PS and SE runs return AJmax to 
about 1 % accuracy is actually quite good. At 
the same time, A J^ax of the FD run is about an 
order of magnitude higher, of about the same level 
as A(J2). 

5.5. Comparison with FD at higher reso- 
lution 

We must point out that in the above compari- 
son, we require the FD runs to use the same grid 
as the SE runs, with similar DOF. This is done in 
order to obtain constrained comparisons that are 
easier to interpret. However, this will, of course, 
make the FD runs intrinsically less accurate, as we 
have so far seen, because FD schemes have a lower 
order truncation than spectral schemes. In reality, 
one can compensate for this and increase the ac- 
curacy of FD runs by using higher resolution or a 
greater number of DOF. This will of course require 
more computational resources. However, with FD 
schemes usually being more efficient {e.g. , scal- 
able), this may certainly be considered a reason- 
able way to obtain more accurate solutions. Here 
we study briefly this possibility by running the FD 
code using higher resolutions. 



In Fig. 17 the blue curve is the fractional er- 
ror AE of the FD run as shown in Fig. I?] (a). The 
green curve is A E' of a FD run using the same grid 
as shown in Fig. [3]but with a 16 x 16 uniform grid 
within each square instead of 8 x 8 as used in the 
blue curve. The purple curve is A E' of a FD run 
using a uniform grid of 256 x 256 resolution. We 
see that the green curve is mostly at a level below 
the blue curve, when the linear resolution is dou- 
bled. The decrease is sometimes quite high, but at 
other times just marginal. For the FD run using a 
256 x 256 uniform grid, which has a cell size equals 
to the smallest cell size of run represented by the 
green curve, the error is substantially (about an 
order of magnitude) lower than that of the blue 
curve. However, when compared with AE of the 
SE run as shown in the red curve in Fig. [t] (a), 
this is still about one to two orders of magnitude 
higher. 

Of course, one can continue increasing the 
DOF of the FD run to try to reach the accu- 
racy level of the SE run. Following the argu- 




Fig. 17. — The blue curve is the fractional error 
AE of the FD run as shown in Fig. |7](a). The 
green curve is AE of a FD run using the same 
grid as shown in Fig. |3] but with a 16 x 16 uniform 
grid within each square instead of 8 x 8 as used in 
the blue curve. The purple curve is AE of a FD 
run using a uniform grid of 256 x 256 resolution. 
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ment in ( [Rosenberg et al. 2006P , the linear error 
bound scaling for the SE case can be written as 
(|Deville et al. 20021 p. 273) 

where p is the polynomial expansion order, h ^ 
1/E is the uniform element length scale, and s is 
the smoothness of the exact solution. For the mo- 
ment we neglect prefactors in the error term. We 
assume s — p so that derivatives can be computed 
up to the order of the method, suggesting a rea- 
sonably smooth function. The error for the FD 
method is 

where N is the total linear (equivalent) number of 
grid cells, and (3 is the nominal truncation order 
(typically k, 2). Equating the logarithmic errors 
yields 

p 

log N = — \og{pE) + prefactor terms. 

This scaling relation provides the number of FD 
cells required to achieve roughly the error of the 
SE method at order p. We cannot actually use this 
scaling directly for comparisons between the FD 
and SE methods, however, because the prefactor 
terms are unknown, and may depend critically on 
the flows. Still, we can illustrate the scaling for the 
case of the SE run in Fig.|7](a). We have that E ~ 
16 (the equivalently-resolved number of elements) 
and p = 8, implying that logA^ ~ 16//3, so with 
/3 = 2, we flnd that achieving comparable accuracy 
in the FD method could lead to the requirement 
for a catastrophically large number of cells, except, 
perhaps, in the case where the filling factor of the 
fine grid is extremely small. This point will require 
further study. 

6. Discussion and Conclusion 

It is worth noting that at the present stage, the 
spectral element method described in this paper is 
more costly in terms of computational time than 
the FD and PS methods with which we compare 
our results, the latter being optimal for periodic 
boundary conditions. We have not tried to com- 
pare CPU usage of the different codes due to the 
fact that the SE and FD codes are running on 
different computational platforms; furthermore, 
both codes are still under continuous development. 



which can change efficiency greatly. Instead, we 
concentrate on comparing accuracy obtained by 
the two schemes, mainly to answer the question of 
whether using SE scheme has any advantage that 
deserves putting more effort into further develop- 
ment. We believe the results of this paper have 
shown that indeed the SE method has substantial 
advantages and needs to be investigated further. 

We have shown that using statically refined 
grids with DOF of roughly linear scaling, SE 
method can produce results with high accuracy 
that can sometimes even approach the spectral ac- 
curacy of the PS method. This can be potentially 
very helpful in the investigation of the 3D MICI 
problem in the ideal limit, which will require adap- 
tive grids that can resolve distinct features. The 
higher accuracy of the SE method can potentially 
provide a reliable definitive answer to the impor- 
tant question of whether a true current sheet forms 
in the Parker's problem of solar and stellar coronal 
heating. 

Our main conclusion that SE methods can pro- 
duce simulations with accuracy somewhat in be- 
tween PS and FD methods is not surprising in 
itself. However, our results yield important quan- 
titative information in the context of statically re- 
fined grids in problems with distinct spatial struc- 
tures. 

The accuracy of FD runs presented here is 
mainly for comparison purposes, since we have 
imposed the restriction of using the same static 
grids with the same DOF. This is by no means a 
suggestion that FD method cannot be used in the 
investigation of the MICI, or similar problems. In- 
deed we have also studied briefly that accuracy in 
FD method can also be increased by using more 
DOF. However, we have not studied the trade off 
of doing this with respect to the increase of the 
usage of computational resources. This may be an 
important topic for future investigation. 
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